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NOTATION 


Unless stated otherwise, the following notation is used. 


a 

A 


e 

pc 

o 


B 

C 

F 

F 

m-n 

h 

i,j,k,l,m,n 

k 

K 

L 

% 

Q 

s 

t 


T 

T c 

f 

V 


Absorptivity 

Area 

Emissivity 
Heat Capacity 

Stefan-Boltzmann Constant (5.66961 (10' 8 ) Wm' 2 K 4 ) 
Time 

Boundary region (B e represents an element boundary.) 
Heat Capacity Matrix 
Force Matrix 

Form Factor between m, n 

Heat Transfer (convective) Coefficient 

Nodes of an element (subscripts) 

Conductivity coefficient 

Conductivity matrix 

Boimdary length 

Flux Normal to Surface 

Internal Heat Generation 

Surface (S e represents an element surface.) 

Thickness 

Temperature 

Ambient Fluid Temperature in convection 

dT/d t 

Volume 


Standard units of the International System have been used in all sample 
problems of this document. 
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FINITE ELEMENT METHOD FOR THERMAL ANALYSIS 

Janeth Heuser 

Goddard Space Flight Center 


INTRODUCTION 

The finite element method (FEM) was applied to structural analysis by Turner et al. about 
17 years ago. As use of the method in the structural field has become widespread and its 
basic ideas are now more clearly defined, the approach is recognized as applicable to any 
problem that can be formulated in variational form. This recognition is significant because 
it stimulates extension of the technique into other fields, such as thermal analysis. Interest 
in this extension has been reinforced by the possibility of coordinating the two areas of 
analysis. Using the FEM approach, structural and thermal data are more easily exchanged; 
using the same computer programs, the two analyses are more closely interrelated and more 
truly representative of things as they are in the real world. 

Conceptually, the FEM involves dividing the continuum into a finite number of elements— 
triangles, for example. The temperature field within each triangle is described in terms of 
the temperatures at the vertices. The heat equation with boundary conditions, written as 
an integral equation in variational form, is minimized. The result is a system of algebraic 
equations which may be solved for the desired solution matrix of predicted temperatures. 

Work has been done to develop and verify the FEM approach for thermal analysis. Compu- 
ter programs exist to handle two-dimensional transient or steady-state heat conduction with 
internal heat generation and the following boundary conditions: convection or a specified 
temperature along the boundary B, and radiation or a constant heat flux normal to B. The 
program described in this paper can handle the two-dimensional case, as well as the 
three-dimensional case subject to the same boundary conditions, or the two-dimensional 
case with radiation, convection, or normal flux over the surface (rather than along the 
boundary) of the two-dimensional continuum. This paper presents the theoretical basis of 
the program, several illustrative problems, and computer program listings. 

THEORETICAL DEVELOPMENT 
Two-dimensional Case 

Consider the equation for two-dimensional (x,y) transient heat flow in a uniformly thick 
plate of thickness t 


3 iv dT \ * 9 ru 9T \ A 3T 

3x x 3x 3y y 3y dr 


0 


( 1 ) 



where the desired temperature T must satisfy the following boundary conditions: 


• Specified temperature T along boundary segment B 

0T 

• Fixed fluxq n normal to B or surface S, q n = 


• Convection between the region and its surroundings along B or over the surface S, 

3T 

Woo -T) = k (— ) 
on 

• Radiation between the region and its surroundings along B or over the surface S, or 
radiation between two boundary or surface segments i and j, oF^fejT* - ocT? ) 


Equation (1) together with its boundary conditions specifies the problem in a unique 
manner. Using the calculus of variations, an alternate formulation is possible by which the 
entire problem may be expressed in one equation— an integral equation in variational form: 


X = 


Iff 


0 s 




3t 


dS dz 


*• t 

If q n T dB dz + JT q n T dS + CC (|hT 2 - liT^T) dB dz 

//(ih T>-KT„T) „S + J'j' q r T dB dz 
S OB 

If * 


0 B 


( 2 )=' 


+ l J q r T dS 
s 


By Euler’s theorem, this function has a minimum represented by Equation (1) with its 
boundary conditions. 

Minimization of X requires an explicit formulation of the temperature T. To find this, 
the region is subdivided into a finite number of triangular elements (Figure 1), and Equa- 
tion (2) may be restated as: X =y]x , where the summation is taken over all elements 
e, and X g represents the functional applied to each triangular element. Assuming the 
temperature field, is hnearly distributed throughout each element, 


within an element. 


T 


a + bx + 


cy = [lx y] 



( 3 ) 


*Note: q f represents radiative flux. The actual definition and treatment of this begins on page 15. 
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Figure 1 . Division of continuum into a finite number of elements. 

To determine the unknowns, a, b, and c, the known temperatures at each of the vertices, or 
nodes, of the triangular element are used to set up a system of equations: 

T. = a + bxj + cy } 

Tj = a + bXj + cyj 
T k = a + bx k + cy k 

In matrix form, this becomes 



Solving the system explicitly by Cramer’s rule yields 

_ ( x jy k - x k yj) T i + ( x k yi - x iy k ) T j + (Vj - x jyi> T k 

1 x i y s 

1 x i y j 

1 x k y k 
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V T j 
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so that 
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where 


2A 


1 

x i 

y i 

1 

x j 

y j 

1 

x k 

y k 


(4) 


is equal to twice the area of the triangle defined by points i ,j, and k. 
Substituting Equation (4) into Equation (3) yields 



This result can be simplified by letting [N] = [1 xyJand[A] = — - times the square 
matrix: ^ 



The element temperature, T, as described by Equation (5) is an explicit and unique linear 
function of the nodal temperatures T j5 Tj, T k and the spatial coordinates x and y. This 
equation provides the desired explicit formulation of T needed to minimize Equation (2). 


3 T 3 T 3 T 

Evaluation of — , — , — with T defined by Equation (5) yields the following explicit 
9x dy or 

expressions for each of these terms: 



( 6 ) 


The functional X=£)X e (Equation 2) can be minimized with respect to the temperature T 
over the entire surface S by the equation 




where Tj , T 2 , • • ■ • T n are the temperatures at the nodal points 1, 2, • • • n in the surface 
S. 
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For each triangular element e, 


dX 

e 

3T 


> = < 


0 

*6* 

5= 

aT. 

0 

• • • 

0 

ax 

e 

3T. 
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0 

0 

ax 

e 

0 
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where i, j. and k are the nodal points of e. 


9X 

Now, consider the expression — defined over the element surface S e with boundary B g . 

8T i ax. ax. 

(Analogous equations will hold for — - and — -.) 

3T j aT k 


ax 

3T 


:■/[//[ 

0 S e 

*f// 


aT_a_3T + ax _a_ ajl 

x ax dT i ax k y 3y STj 3y 
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aT A* A l 

pc — — dS dz + 
at 3T : 
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e 
t 


ff 


q dB dz + 
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ff 

S e 
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Q — dS 
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Vf 

0 B 

e 

ff ^fr dB it * ff q, fjr 


ff 
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3T 
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aT 
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( 7 ) 


d X 

To evaluate — ^ , consider each portion of the integral equation separately. First, consider 
dT i 

the double integral over S„ dealing with conduction. 
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Conduction 


0 X 0 X 

Note that by defining — and — explicitly as in Equation (6) 
0x 3y 


d = [0 10] [A] { 0 ) = a 21 = a.j 


9Tj 9x 


3 3T 


3Tj 3y 


^ = [0 0 1] [A] 0 = a 


31 a ki 


Thus 

t 


///[• 

0 S e 


3T 3 3T 


+ k. 


3T 3 3T 


-fff 


x 3x 3T . 3x y 3y 3T; 3y J 

,1 


dS dz- 


X 


[0 1 0] [A]/ Tj 



a ji + k y 


[0 0 1] [A] 



a ki 


= fff [ Rx (3jiTi + ajjTj + 3 i k T k )a ji + k y ( a ki T i + a kj T j + a kk T k ) 


Assuming constant k x and k y over the element 


-fff 


t a ji a ki] 


" l a ji a kil 


k x 0 

o k. 


a ji a ii a .ik 
a ki a kj a kk 


k x 0 


0 k. 


a., a-- a-, 
Ji JJ Jk 


a ki a kj a kk 



Tj ) dx dy dz 


dx dy / dz 


/' 


= d^yj-y k x k - x jJ 


r k 

0 

1 

X 

0 

k y _ 

2 A 


yj - y k yk - y s yi - ^ 

x k- x j x i - x k x j - x i 


dx dy dz 


i ki dx dy dz 
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/ 3 X CONDl t r 

' k x 

0 - 

' y j - y k 

y k “ y i 

m 

M 

( 9Tj r 4a y i x il 

_ 0 

k y- 

- x k- x j 

x i - x k 


T i 

1 T J 


Combining this with similar expressions for ^ ON - and ^° ND results in 

3T i 3T k 



( 8 ) 


Performing these operations for each element and summing the resulting matrices yields the 
total n-by-n conductivity matrix K for the continuum (\yhere n = the number of individual 
nodes). The above portion of the K matrix will contain nonzero terms to be added to K u , Ky, 
K ik ,K.., K jj> K jk , K ki , K kj) K kk , and zeros elsewhere. ^ 


Internal Heat Generation 

The internal heat generation portion of this equation is 


dX, 


0T- 


dS dz = 


■ -///<% 

* 0 S e 0 

= jy Q ( a u + a ji x + a kiy) dx dy J dz 


fff Q [i xyj [A] jo 


dx dy dz 


Assuming constant Q over the element 
= -Q A (a ;i + aji x + a kj y) t 


since 


If 


dx dy = A, x = ^ x dx dy, y - y dx dy 


where x and y are coordinates of the centroid of the element. Now, x =(l/3)(Xj + Xj + x k ) 
and y =(l/3)(y j + y- + y k ). Substituting these into the equation yields 
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L ii K I 

a ii + y( x i + x j + x k ) + T" (y i + y i + y *> 

QAt ( (y j _y k> 

= I x. y, - x. y. + - - (x. + x. + x, ) 

2A \ J k Uc J 3 v i J 


O k - Xj) 


( y i + Vj + y k ) 


1 x i y i 

Q At , 1 . , 

— ( 2A> 1 X i y j 

1 X k y k 


9^1 ( J_ ) (2A) = - 9^1 
3 1 2A J K ' 3 


ew , . 3X Q _ aX Q _ QAt 
Similarly, ^ ^ r 

the force matrix is given by 


Thus if F =y> for each element, the contribution to 


— Q A t/3 

— Q A t/3 

— Q A t/3 


where -QAt/3 is added to the terms F ; , F ; , and F. . 

: 1 J K 

Heat Capacity j j 

The contribution due to the heat capacity term is 


-j = JJJ pc££-dSdz = JJJ pc [1 x y] [A] \ Tj \ [lxy] [A]lo\ dx dy dz 

1 'e OS 1 n / t \ / 0 \ 


Assuming constant pc over e 


= P c If ^ a ji x + a ki y)[l X y] [A] < Tj > dxdyj dz 

W 0 

= P c JJ [ a ii + a ji x + a ki y a ii x + a ji x2 + a ki xy a ii y + a ji xy + a ki y2 


dx dy [A] /Tj \ t 
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Using the identities £j'* 2 dx d y = ^( x i x j + x i x k + x j x k + E x m ) and JJ ' x y dx dy 

= 4(E x m y n + E x m y m )’ where m = h u k and n = h u k ’ this becomes 


3X 


3T : 


= pet 


+ 2a.. 

m ji 


A f 1 2a ii + 4 aji^ x m + 4 a ki 5^y m , 4 * 0^2 

1 z L m m m 

( x i x j + x i x k + x j x k + a kiE X m y n +a k iE X m y ^ 4a iiS y 


m 


+ a ; 


i]C x m y n + 

a ji 4 

m,n 


a. + a. x 

+ a 

im jm n 



Defining bj^ = a im +a j m x n + a km y n , this may be put into the form 


f 3X 


3T 


pc\ = pet 

7T 


"I 

'2 

1 

r 


b- b( bfi 

1 

2 

i 


L J 

_1 

1 

2 _ 



b l bj b> 

b i b ] bJ k 
bf bf b k k 


{t}= eg. [2 , r 


3X 3X pc 

for the two-dimensional triangle. Consideration of — — and -r= — shows that the 

9T j aT k 

additions to the total heat capacity matrix [C] due to this element are 

f a V\ _ pctA 

1 TT L TT 


2 

1 

1 


T i) 

1 

2 

1 

< 

T j 

1 

1 

2 




(10) 


Thus far, three matrices (K, F, and C) have been defined. The conductivity matrix, K, is 
multiplied by the vector 



the force matrix, F, is multiplied by one; and the heat capacity matrix, C, is multiplied 
by 
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In the evaluation of the remaining integrals of X, if a term is to be multiplied by | T it 
will be placed in the appropriate row and column of [K ] . If it is neither a coefficient 
of JT} nor of | T}, it will be placed in { F } . 


Convection 

Case 1 : Convection along boundary of thickness t — For convection along an element 
boundary (i, j), one again assumes that T varies linearly so that the temperature T at any 
point on the boundary of length L is given by 

T = Tj + (Tj - T.) £ 

Now 


3T\ 


//( 


hT — - hT-o | dB dz 


3T- 


3T \ 
3Ti ) 


0 B e 
L 


/ *>[t, ♦ (Tj ~ Tj) l] [(i - £) - d x y d 2 


0 

If h and are constant over B„ 


= h 


(Tj - 2Tj + T 00 )L 2 (Tj - T.)L* 

(Tj - Too)L + — + J 


2L 


3L 2 


T T T. T- 

, 1 1 00 *i 1 

= hLt l T. - + — - T ; + - — - + - 

1 00 2 1 2 3 3 


- hL '(j T ‘ + * T J " 5 T -) 


= [ hj f ^ °] - ' hLT ” 1 ' 2 


0 


dx h ax h 

The expression for is similar, and = 0 because T k is not a part of the convec- 

j k 


tive boundary. Thus 


1 3T r 


hLt/3 hLt/6 


LhLt/6 hLt/3 J LTjJ ^hLT^t/2 


‘F 

i J v- 


hLT^t/2 


(ID 


where the terms of the square matrix are added into [K], increasing K ii; K^, K^, and Kjp 
and the terms of the column matrix are added to F- and F.. 



Case 2: Convection over surface — Now assume that convection exists not only along an 
outside boundary but over the entire element surface (ij, k). Then 



Assuming h and constant over the element surface 



hA 

12 


[2 1 1 ] 



hA 

3 


T 


OO 


The terms 
of 



and 


9X h 

— - are found in a similar manner, yielding a final convective matrix 
3T k 



2 

1 

1 



fhAT^/3] 
< hAT^/3) 
^hAT^/3] 


( 12 ) 


Each element of the 3-by-3 matrix (each term multiplied by hA/12) is added to its corres- 
ponding (i, j, k) component of the larger matrix [K] . Likewise, — hTooA/3 is added to Fj, Fj, 
andF k . 
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Specified Flux Normal to Boundary or Surface 

Case 1: Normal to boundary of thickness t — In this situation 


( _li]= 1/ a dBdi 

\ 3 T J JJ " 9T i 

\ /e 0 B e 


Since T is assumed to be linear over the element, it also is linear along a boundary (i, j) and 
thus may be written as T = Tj + (Tj - Tj) x/L, where L is the length of (i, j). Using this 


expression 


dT _ x 
3T. L 


3T 

3T~ 


x 

L 


3T 

ST,. 


= 0 


Therefore, assuming that the flux is normal to boundary (i, j) 
L 

r 

dz 


(— f q (1 - -) dx f 

\ 3T. / J n V L 7 J 
\ */e 0 0 


• f q n is constant on (i, j) 


‘ % ft) 


L 

• - <!.«•-?>• - — 


Similarly 


— ' « f q * dx f dx = (q *—)l 

3 Tj J 4 " L J 2 L ; 


q n Lt 


3X 


3T. 


q = 0 
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The flux matrix for this element is, therefore 



tq n L/2*) <- i 
tq„L/2 ? «- j 
0 J «- k 


(13) 


Since this matrix is multiplied by neither 
terms Fj and Fj by (tq n L/2). 


{t} nor {i} 


it is added to 


w 


, increasing the 


Case 2: Normal to surface — If the flux is normal to the surface of element (i, j, k), then 



where q n is constant over (i, j, k) 

= % ff 3ii + 3j ‘ X + akiV ^ 

q n A 

= — ( 3a ii + a ji ( x i + x j + x k) + a ki + yj + y k » 




< 2a > - 


<TA 

3 


Likewise, 8 X q /3Tj = 3X q /9T k = q n A/3, and the flux addition to the |f| matrix becomes 


{ 



f<ln A / 3 ) 

q n A /3 

Un A / 3 J 


(14) 


where q n A/3 is added to F., F., and F k . 
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Radiation 


Radiation between element e and other radiative elements R or between boundary e and 
other radiative boundaries R may be described for the transient situation by 

q, - £ (oe e A,fVR - »<*eMe-R 

R 

where e, a, and F „ are assumed constant over individual elements 

5 5 e-R 

= ae e A e T e XX-R ffa e A eE F e -R T R 


= M e A e T !^ fl e A «J F e-R T R 

~ <Je e A e 1 (e)- <Ta e A eXI F e-R 1 (|l 


where T^is the temperature of boundary e found in the previous time step of the 
computer-iterative algorithm. This approximation is numerically stable and accurate if 
5 r is chosen small enough to ensure small values of | T 4 — lA | and I - Tj^ |. Treating 
T 4 in this manner leads to the explicit evaluation of a q r at tne beginning of each time 
step. Hence, during that time step, q r may be treated as a specified constant flux: 


(x,) = ff (« e l^-o« e X;F e . R T^)IdS 
e sf 



S e 


(15) 


Case 1: Radiation between two boundaries or between a boundary and a constant 
temperature heat source — Letting q n = (ne e T 4 - oa e yi F e-R T ^’ this case is handled in 


R 

the same way as Case 1 of “Flux Normal to Boundary or Surface,” the previous section. 
The expression corresponding to Equation (13) is 



!(«eT©-™«Z F e-R T ® l > i 

0 J <- k 


(16) 
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where T 0 = (^+ 2 when i and j are nodes of the radiating boundary e. Also, =(T 0 

+ T @)/2 when is on an element boundary, (p, m) and = T R when T R is a specified- 
temperature heat source. 

Case 2: Radiation between two surfaces or between a surface and a constant-temperature 
heat source — Following the procedure in Case 2 of the previous section, the following 
addition to the F matrix results 


^ V 3 (° e e 1 (t)- oa eZI F e-R 1 |))h i 

V 3 < 17 > 

V 3 


where (^+ T 0)/3 and i, j, and k are nodes of the radiating element e. Also, 

T 0 = (T^p + T^gj + T 0>/3 if T r is a part of the modeled continuum, and = T R if T R 
is a specified-temperature heat source. 


Three-dimensional Case 

The previous development for two-dimensional heat transfer may be expanded to handle 
the three-dimensional heat equation 


_ 3 _ 

dx 


< k x 


!I) + i- (k *1 

dx J dy y y dy 


) 



9T 

3z 


+ 


Q - pc 


9T 

dr 


0 


(18) 


where the desired T must satisfy the same boundary conditions given for the two-dimensional 
case. 

In variational form this problem becomes 


X = 


m 

V 

'// 




dV + 


fU 


T dS 


hT - hT TO T dS + 


/A 


(19) 


TdS 


which, as before, may be minimized to yield |T j . Assuming that T is linear over the 
three-dimensional tetrahedron and proceeding as for two dimensions, one may derive 
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T = a + bx + cy + dz = [1 x y z] 




= [ix yzl j- [A] {t} 


for the tetrahedron 


where 

1 x i *i z i 

1 x i y j z j 

1 x k y k z k 

1 x i yi z i 

V = t is the volume of the tetrahedron. 



x j 

y j 

z j 


1 

y j 

z i 

a i = 

x k 

y k 

z k 

b i = 

1 

y k 

z k 


X 1 

*1 

Z 1 


1 

yi 

z i 



x j 

i 

z j 


x i 

y j 

i 

c i = 

x k 

i 

z k 

d i = 

x k 

y k 

i 


x i 

i 

z i 


x i 

yi 

i 


(20) 


( 21 ) 


Other constants may be defined by cyclic interchange of the subscripts in the order ij,k,l. 
Each time a cyclic interchange is made, it is necessary to multiply the newly defined constant 
by -1. 
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The process of minimizing X e with respect to T leads to 


3X 

3T. 


, ax a aT , aT a 3T 

k — - — — — + k. 


Iff\ x ax aTj 3x ' "y ay dT. a y J 

- Iff » Sr " • Iff "" 


dV 


at 3T; 


dv + ff 0 . If ds (22) 


ff 


t ,J (KT 

s, 


- hT M dS + 
aT. 00 3T. 


ff aT 

JJ ^ 3T. 


dS 


Note that the radiative, flux, and convective terms in this equation are identical to those 
for the surface conditions in the two-dimensional analysis. Hence, their evaluation has been 
completed. Now, consider the evaluation of the conductive, internal heat generation, and 
heat capacity integrals. 


Conduction 

Following the same procedure as for two-dimensional conduction, but including k z and 
redefining T, the addition to the [K] matrix for element (i, j, k, 1) becomes 



Nonzero terms are added to [K ] where m = ij,k,l and n = ij,k,l. 

Internal Heat Generation 
Using 



and 


( 22 ) 


y 

z 


1 _ 

V 

\_ 

V 
1 _ 

V 


/ 

/ 

/ 


x dV (tetrahedron) 
y dV 
■l dV 
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where (x , y , z ) are coordinates of the centroid of the tetrahedron, and following the method 
used for internal heat generation in two dimensions yields 



-QV 

where — - is added to F . , F . , F. , and F. . 

4 1 J K i 

Heat Capacity 

The contribution due to the heat capacity term is evaluated by using the identities 


f x2 dv = To \ ( E } + x i x j + x i x k + x i x i + x j x k + x j x i + x k x i 1 

J L m J 

/* XydV = IU tE X m y n + E X mym] 

where m = i, j, k, 1; n = i, j, k, 1. 

Again, by letting b” = (a m + b m x n + c m y n + d m z n ), the three-dimensional equivalent 
of Equation (10) becomes 


m- 


/pcV 

V 20 


b} b] bf b? 
bt bi b^ bj 


J 


i J 


J 




b i b i b k 
b i b i »i 


2 111 
12 11 
112 1 
1112 


2 111 
12 11 
112 1 
1112 


b i b ; 


i J 

bf bf 

b! b J 
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General Matrix Equation 


Performing the above operations for each element and summing the resulting matrices results 
in a system of simultaneous linear equations which may be written in matrix form as 

[K]{t} + [C]{t} + {f} - 0. (25) 

To facilitate computer solution, time derivatives are replaced by finite difference approx- 
imations, yielding 


( I k t 1 + f- [C]){t t } - (j-l C] - [K t _ 5t )){t t _ 5 ,} - {f 7 + F,_ 5t } 

From Equation (25) 

[K r ]{T r } t [C J {t ♦ {f,} = 0 

[K t _ St ){t t _ 5t } + [C, {T r _ Sr } + {F r _ Sj }= 0 


Then 




or 


N = l <{ T r}-{v4 -{b-5,} 

Substituting this result into Equation (25a) gives 

[ ^'H + [C1( MTr}-{T T _4>- {b- 5 r}) * { F r} = ° 


or 


([K rl + b I C )){t t } = [C] ({t r _ Sr } + i- ) -{f 7 } 


(26) 

(25a) 

(25b) 


( 27 ) 


Now, from Equation (25b) 

{Vi,}= l C r‘ (-l K r-Jrl{T r - 6 r }-{ F r-Sr}) 
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Then, substituting this into Equation (27) yields 


«M + ^rfd) { T r} - [CK-tC]' 1 [K r _j T ]{T T _j r }- [C]->{f t _ 6t }) 

* [C] h { T r- 5 r}-{^} 

which can be reduced to Equation (26). To find only the steady state distribution for a 
problem with no radiative boundaries, the matrix equation [K] W + F} = 0 is used. 
Solution of this linear system yields the steady state results in one step. Using Equation (26) 
instead would require a number of 5 r iterations before reaching the steady state tempera- 
tures. For problems including radiation, it is necessary to use Equation (26) with an approp- 
riate 5r and iterate to the steady state temperature distribution. (This iteration is necessitated 
by the use of the approximation T* = T^.) 

A computer program has been written to solve the matrix Equation (26). The program also 
includes consideration of the boundary condition that holds specified boundary nodes at 
constant temperatures. To consider such a boundary condition, let 

A - [K r ] ♦ £[C] 

B - £tC]-[K T _ Sf ] 

C - -<{ F r} + K-4>' 

Then 

A M ■ b{ T t- 5r} + C. (28) 


Suppose that of the n nodes representing the thermal model, the last k nodes are held fixed. 
Thus the first m = n — k nodes are variable. The matrix Equation (28) can be partitioned 
into m-by-m and k-by-k submatrices as follows. 


L m b y m 


k by m 


l m by k 


l k by k 



® m b y m b y k| 


k by m 


B 


k b y k 


tWr-Sr 


(29) 
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Matrix multiplication yields two matrix equations 

+ ^mk^k “ B mrn^(m) + ®mk^k) + (29a) 

A km T m + A kk T k = B km T © + B kk T © + C k < 29b ) 

Only the first equation need be solved for the unknown temperatures |T m |. If the specified 
nodes are not necessarily the last k, a similarity transformation may be applied to cast the 
system into this form. After solution is obtained, again it may be applied to return the 
answers to their original order. 

Sample Problems 

Following are some sample problems in one, two, and three dimensions to which the pre- 
ceding method was applied. Each sample includes a description of the problem and the 
temperature distributions predicted at certain times as well as at steady state. For the 
first three problems, exact or finite difference temperature predictions also are given for 
comparison. 

Problem 1: One Dimension— Finite Element and Exact Solutions 

One-dimensional convection occurs at the surface of a semi-infinite solid with constant 
initial temperature T Q into a medium at 0°K. The exact solution is given by 

= erf — + e hx + h 2 ar er fc [— + harl 
Tq 2 ar l2ar J 

where a = k/pc is the thermal diffusivity. 

Letting k = 1 Wm^K' 1 , pc = 1 J m" 3 K _1 h = 1 Wm^K' 1 , T Q = 100°K and = 0°K, the 
temperature distributions (exact and finite element) for various times are plotted in 
Figure 2. As can be seen, agreement of results is excellent. 1 


'Problem and analytical solution graphs were obtained from the work by R.V.S. Yalamanchili and S.C. Chu 
(see Sources). 
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The following thermal properties are used in sample problems 2 to 6. 

Material: aluminum 

k x = k y = k z = 62.99 w m l K" 1 

pc = 72.76 Jm^r 1 

- h = 42.03 W m‘ 3 K -1 ' 

e = 0.1 = a; absorptivity 

T 0 = 297.2K 

t = 0.0254 m 

Problem 2: Two Dimensions — With Radiation and Convection (Finite Element and 
Finite Difference Solutions) 

Consider one quarter of a symmetric pipe where the external surface, a, radiates to an 
environment at 255. 6K and the internal surface, b, is a convective boundary with 
= 373.3K. 


1.2192 m 



20 . 21 22 



255.6 K 


Table 1 shows transient and steady state temperature predictions for nodes 1 to 23, found 
by the finite element program. Also shown are ^predictions for corresponding nodes of an 
implicit finite difference program using a 30-node model. 
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Table 1 


Comparison of Finite Element and Finite Difference 
Results for Problem 2. 


Node 

T 

at 1800s 

T 

at 3600s 

T 

at 7200s 


Finite 

Element 

Finite 

Difference 

Finite 

Element 

Finite 

Difference 

Finite 

Element 

Finite 

Difference 

1 

297.4K 

297.9K 

300.6K 

301. OK 

308.6K 

308.7K 

2 

298.0 

298.4 

301.7 

301.9 

309.7 

309.7 

3 

298.6 

298.9 

302.7 

302.9 

310.8 

310.7 

4 

299.0 

299.1 

303.3 

302.3 

31 1.3 

311.1 

5 

297.6 

298.1 

300.9 

301.3 

308.9 

308.9 

6 

298.1 

298.4 

301.6 

301.8 

309.6 

309.6 

7 

299.1 

299.3 

303.0 

303.2 

310.9 

310.9 

8 

299.9 

300.2 

304.3 

304.4 

312.2 

312.2 

9 

300.2 

300.4 

304.8 

304.9 

312.7 

3-12.7 

10 

298.8 

299.2 

303.1 

303.1 

311.0 

310.9 

11 

299.9 

300.1 

304.2 

304.2 

312.1 

311.9 

12 

. 302.5 

302.5 

306.9 

306.9 

314.6 

314.4 

13 

304.4 

304.5 

309.1 

309.2 

316.6 

316.6 

14 

304.9 

305.0 

309.8 

309.9 

317.2 

317:3 

15 

300.7 

300.7 

305.6 

305.3 

313,4 

313.1 

16 

302.7 

302.3 

307.5 

307.2 

315.2 

314.8 

17 

309.6 

309.2 

313.9 

313.6 

320.9 

320.1 

18 

313.4 

313.6 

317.8 

317.8 

324.3 

324.4 

19 

314.1 

314.4 

318.6 

318.9 

325.1 

325.3 

20 

301.7 

301.6 

306.8 

306.6 

314.6 

314.3 

21 

303.9 

304.1 

308.9 

308.9 

316.6 

316.4 

22 

311.9 

312.4 

316.5 

316.8 

323.2 

323.4 

23 

297.1 

. 297.7 

•299.9 

300.6 

307.9 

308.2 
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Table 1 (Continued) 


Comparison of Finite Element and Finite Difference 
Results for Problem 2. 


Node 

T at 18,000s 

T at 36,000s 

T for Steady State 

Finite 

Element 

Finite 

Difference 

Finite 

Element 

Finite 

Difference 

Finite 

Element 

Finite 

Difference 

1 

328.1K 

327.8K 

347.9K 

347.7K 

370.9K 

370.9K 

2 

328.9 

328.7 

348.4 

348.7 

371.1 

371.1 

3 

329.1 

329.4 

348.8 

348.6 

371.1 

371.1 

4 

330.1 

329.7 

349.1 

348.7 

371.2 

371.2 

5 

328.3 

328.1 

348.1 

347.8 

370.9 

370.9 

6 

328.8 

328.5 

348.4 

348.1 

371.1 

371.1 

7 

329.8 

329.5 

348.9 

348.7 

371.2 

371.2 

8 

330.7 

330.4 

349.4 

349.2 

371.3 

371.3 

9 

331.0 

330.8 

349.7 

349.4 

371.3 

371.3 

10 

329.8 

329.4 

348.9 

348.6 

371.1 

371.1 

11 

330.6 

330.2 

349.4 

349.1 

371.2 

371.2 

12 

332.3 

331.9 

350.4 

350.1 

371.4 

371.4 

13 

333.7 

333.6 

351.2 

351.0 

371.4 

371.5 

14 

334.2 

334.1 

351.4 

351.3 

371.5 

371.5 

15 

331.6 

331.0 

349.9 

349.5 

371.2 

371.2 

16 

332.8 

332.2 

350.7 

350.2 

371.3 

371.3 

17 

336.8 

336.1 

352.9 

352.4 

371.6 

371.6 

18 

339.2 

339.1 

354.2 

354.1 

371.7 

371.8 

19 

339.7 

339.6 ' 

354.6 

354.4 

371.8 

371.8 

20 

332.3 

331.9 

350.3 

350.0 

371.2 

371.2 

21 

333.7 

333.4 

351.2 

350.9 

371.4 

371.4 

22 

338.4 

338.3 

353.8 

353.7 

371.7 

371.7 

23 

327.6 

327.5 

347.6 

347.4 

370.8 

370.8 
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Problem 3: Two Dimensions, With Internal Heat Generation, Specified Flux, and 
Convection (Finite Element and Finite Difference Solutions) 

Consider one quarter of a symmetric pipe with the external boundary, a, exposed to the 
normal flux described by the curve below; with the internal convective boundary, b, 
exposed to a fluid at 255.6K; and with internal heat generation throughout the pipe, 

(Q = 10.34 W m' 3 ). The model used is the same as that for problem 2, illustrated previously. 
Results are shown in Table 2. 


1.2192 m 




1 

0.6096 m 

Q = 10.34 


255.6K 



| 1.2192 m 


Flux for boundary a 



T,s 


27 



Table 2 


Comparison of Finite Element and Finite Difference 
Results for Problem 3. 


Node 

T 

at 1800s 

T 

at 3600s 

T 

at 7200s 

Finite 

Element 

Finite 

Difference 

Finite 

Element 

Finite 

Difference 

Finite 

Element 

Finite 

Difference 

1 

298. 8K 

298.2K 

298. OK 

297.5K 

295.2K 

294. 8K 

2 

298.2 

297.8 

297.0 

296.7 

294.1 

293.9 

3 

297.7 

297.4 

296.2 

296.0 

293.3 

293.1 

4 

297.5 

297.3 

295.8 

295.8 

292.9 

292.8 

5 

298.7 

298.2 

297.8 

297.3 

295.0 

294.6 

6 

297.8 

297.4 

296.8 

296.5 

294.0 

293.7 

7 

296.9 

296.7 

295.6 

295.4 

292.8 

292.6 

8 

296.3 

296.2 

294.7 

294.6 

291.8 

291.6 

9 

296.1 

296.0 

294.3 

294.2 

291.4 

291.3 

10 

297.7 

297.4 

296.3 

296.1 

293.4 

293.3 

11 

296.5 

296.3 

295.0 

294.9 

292.2 

292.1 

12 

294.7 

294.7 

293.1 

293.1 

290.3 

290.2 

13 

293.6 

293.5 

291.6 

291.5 

288.9 

288.7 

14 

293.2 

293.2 

291.2 

291.1 

288.4 

288.2 

15 

296.6 

296.5 

294.7 

294.8 

291.8 

291.9 

16 

294.8 

294.9 

292.9 

293.1 

290.1 

290.3 

17 

290.8 

291.1 

289.0 

289.3 

286.4 

286.7 

18 

288.5 

288.4 

286.7 

286.4 

284.2 

284.1 

19 

288.1 

288.0 

286.1 

285.9 

283.7 

283.5 

20 

296.0 

295.9 

294.0 

294.0 

291.1 

291.1 

21 

294.1 

294.1 

292.1 

292.1 

289.3 

289.3 

22 

289.4 

289.2 

287.6 

287.4 

285.1 

284.9 

23 

299.6 

298.6 

299.0 

298.0 

296.2 

295.3 
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Table 2 (Continued) 


Comparison of Finite Element and Finite Difference 
Results for Problem 3. 


Node 

T at 

18,000s 

T 

at 36,000s 

T for Steady State 

Finite 

Element 

Finite 

Difference 

Finite 

Element 

Finite 

Difference 

Finite 

Element 

Finite 

Difference 

1 

288. IK 

287. 7K 

280. 7K 

280.3K 

272. 2K 

272.4K 

2 

287.1 

286.9 

279.9 

279.6 . 

271.4 

271.8 

3 

286.3 

286.2 

279.2 

279.0 

271.0* 

■ 271.3 

4 

286.1 

285.9 

279.0 

278.8 

270.8 

271.2 

5 

287.9 

287.6 

280.6 

280.2 • 

272.1 

272.4 

6 

286.9 

286.7 

279.7 

279.4 

271.3 • 

271.7 

7 

285.9 

285.7 

278.8 

278.6 

270.6 

270.9 

8 

285.1 

284.9 

278.1 

277.9 

270.0 

270.3 

9 

284.7 

284.6 

277.8 

277.6 

269.8 

270.1 

10 

286.5 

286.4 

279.4 

279.2 

271.2 

271.6 

11 

285.4 

285.3 

278.4 

278.2 

270.3 

. 270.7 

12 

283.8 

283.8 

277.1 

277.2 

269.3 

269.7 

13 

282.6 

282.4 

276.1 

275.9 

268.6 

268.9 

14 

282.2 

282.1 

275.8 

275.6 

268.4 . 

268.7 

15 

285.2 

285.2 

278.4 

278.3 

270.5 

270.3 

16 

283.7 

283.8 

277.1 

276.9 

269.4 

269.8 

17 

280.7 

280.9 

274.7 

274.7 

267.8 . 

268.2 

18 

278.8 

278.7 

273.2 

273.1 

266.7 

266.9 

19 

278.3 , 

278.2 

272.8 

272.7 ■ 

266.4 

266.7 

20 

284.6 • 

284.6 

277.9 

277.8 

270.2 

270.6 

21 . 

283.1 

283.0 

276.6 

276.4 

274.7 

269.4 

22 

279.6 . 

■ 279.6 

273.9 

273.7 

267.2 

267.3 

23 

289.0 

288.2 

281.6 

280.7 

272.9 

273.3 
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Problem 4: Three Dimensions - Specified Temperature and Convective Surfaces 

Consider a hollow cube where the outside surface temperature for side a is held constant at 
422.2K and the convective surface of side b is exposed to a fluid with = 255. 6K. Tetra- 
hedral elements have been used to model an eighth of this symmetric cube, resulting in the 
nodal points shown in the sketch. Results of the thermal analysis are shown in Table 3. 


SURFACE THICKNESS 0.0254 m 




POINTS ON INSIDE SURFACE 


POINTS ON OUTSIDE SURFACE 




Table 3 

Finite Element Results for Problem 4. 


Node 

T at 120s 

T at 300s 

T at 900s 

T at 1800s 

T at 3600s 

T for 
Steady State 

1 

422. IK 

422.2K 

422. 2K 

422.2K 

422.2K 

422.2K 

2 

422.1 

422.1 

422.2 

422.2 

422.2 

422.2 

3 

416.9 

418.5 

419.7 

420.1 

420.3 

420.3 

4 

413.9 

416.5 

418.4 

419.1 

419.3 

419.4 

5 

366.2 

383.9 

396.4 

400.8 

402.9 

403.3 

6 

323.8 

350.3 

371.1 

379.2 

383.1 

383.8 

7 

321.8 

348.1 

369.9 

378.4 

382.4 

383.2 

8 

303.1 

323.7 

346.3 

357.1 

362.3 

363.3 

9 

296.2 

305.2 

323.1 

334.6 

340.6 

341.7 

10 

296.5 

305.8 

325.2 

337.2 

343.2 

344.4 

11 

294.4 

293.9 

305.2 

316.3 

322.2 

323.3 

12 

290.3 

284.8 

289.3 

298.3 

303.3 

304.2 

13 

292.4 

288.1 

295.6 

305.9 

31 1.5 

312.6 

14 

287.6 

281.1 

283.1 

290.8 

295.3 

296.1 

15 

285.6 

277.3 

275.6 

281.7 

285.4 

286.2 

16 

422.2 

422.2 

422.2 

422.2 

422.2 

422.2 

17 

422.2 

422.2 

422.2 

422.2 

422.2 

422.2 

18 

422.2 

422.2 

422.2 

422.2 ’ 

422.2 

422.2 

19 

422.2 

422.2 

422.2 

422.2 

422.2 

422.2 

20 

365.9 

383.8 

396.3 

400.7 

402.8 

403.2 

21 

323.6 

350.3 

371.1 

379.3 

383.1 

328.3 

22 

322.1 

348.1 

, 369.7 

378.2 

382.2 

383.0 

23 

303.1 

323.7 

346.3 

357.1 

362.3 

363.3 

24 

296.2 

305.3 

323.1 

334.7 

340.6 

341.7 

25 

296.7 

306.0 

325.5 

337.5 

288.0 

344.7 

26 

294.4 

294.0 

305.4 

316.4 

322.3 

323.4 

27 

290.1 

284.2 

288.3 

297.2 

302.2 

303.2 

28 

292.3 

287.4 

294.6 

304.9 

310.6 

31 1.6 

29 

287.5 

281.0 

282.9 

290.7 

295.1 

295.9. 

30 

285.2 

277.1 

275.2 

281.3 

285.0 

. .285.7 
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Problem 5: Three Dimensions - Specified Temperature, Flux, and Radiative Surfaces 

The outside surface temperature for side a of a hollow cube is specified as T a = 233. 3K; 
the shaded portion of the cube is exposed to normal flux, an external q = 1398.1 W m' 2 ; 
and radiation exists between all internal surfaces. Tetrahedral elements have been used to 
model the entire cube. The grid point locations of the models are shown in the “flattened 
diagrams.” 
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Results of the finite element analysis are shown in Table 4. The temperatures for nodes 
27,30, 33, 36, 39, 42, 45, 48, and 51 are constant (T = 233. 3K) and hence are not shown 
in this table. 
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Table 4 


Finite Element Results for Problem 5. 


Node 

T at 1 80s 

T at 300s 

T at 900s 

T at 3600s 

T at 7200s 

T for 
Steady State 

1 

236. IK 

235. IK 

234.5K 

233.9K 

233.6K 

233.4K 

2 

286.6 

266.6 

256.4 

245.3 

237.5 

235.3 

3 

298.3 

284.5 

269.4 

252.2 

239.9 

236.4 

4 

235.6 

234.8 

234.3 

233.8 

233.5 

233.4 

5 

285.9 

266.7 

256.7 

245.6 

237.7 

235.4 

6 

297.6 

285.5 

270.4 

252.8 

240.3 

236.8 

7 

236.7 

235.4 

234.8 

234.1 

233.6 

233.4 

8 

286.7 

266.9 

257.0 

245.8 

237.9 

235.7 

9 

298.1 

284.6 

269.9 

252.7 

240.4 

236.9 

10 

235.6 

234.8 

234.3 

233.9 

233.6 

233.4 

11 

285.2 

267.1 

257.6 

246.6 

238.8 

236.6 

12 

297.5 

286.7 

272.2 

254.7 

242.3 

238.7 

13 

235.9 

235.0 

234.6 

234.0 

233.7 

233.6 

14 

287.2 

269.3 

259.9 

249.0 

241.2 

239.0 

15 

300.5 

289.4 

275.2 

258.1 

245.8 

242.3 

16 

235.6 

234.9 

234.4 

234.0 

233.7 

233.6 

17 

287.3 

269.4 

259.9 

248.9 

241.2 

239.0 

18 

298.4 

288.2 

273.8 

256.4 

244.0 

240.5 

19 

236.7 

235.5 . 

234.9 

234.2 

233.8 

233.6 

20 

286.6 

267.4 

257.6 

246.6 

238.8 

236.5 

21 

297.9 

285.3 

270.8 

253.8 

241.7 

238.3 

22 

235.6 

234.8 

234.3 

233.9 

233.4 

233.4 

23 

285.8 

266.5 

256.7 

245.7 

237.9 

235.7 

24 

297.4 

285.4 

270.5 

252.9 

240.6 

237.0 

25 

233.6 

233.5 

233.4 

233.4 

233.3 

233.3 

26 

295.8 

231.8 

272.4 

254.3 

241.6 

237.9 

28 

286.6 

266.6 

256.4 

245.3 

237.5 

235.3 

29 

299.1 

285.5 

270.2 

252.5 

240.1 

236.5 

31 

285.9 

266.8 

256.7 

245.6 

237.7 

235.4 

32 

298.1 

286.2 

270.9 

253.1 

240.4 

236.8 

34 

286.8 

266.9 

257.0 

245.8 

237.9 

235.7 

35 

298.7 

285.3 

270.4 

252.9 

240.5 

236.9 

37 

285.3 

267.2 

257.6 

246.7 

238.8 

236.6 

38 

298.1 

287.4 

272.8 

255.1 

242.4 

238.9 

40 

287.3 

269.4 

260.0 

249.1 

241.4 

239.2 

41 

301.5 

290.8 

276.3 

258.8 

246.4 

242.8 
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Table 4 (Continued) 

Finite Element Results for Problem 5. 


Node 

T at 180s 

T at 300s 

T at 900s 

T at 3600s 

T at 7200s 

T for 
Steady State 

43 

287.3K 

269.4K 

259.9K 

249. OK 

241.2K 

239. OK 

44 

298.9 

288.9 

274.3 

256.8 

244.3 

240.7 

46 

286.7 

267.3 

257.5 

246.5 

238.7 

236.4 

47 

298.6 

286.1 

271.4 

254.1 

241.8 

238.3 

49 

285.9 

266.6 

256.7 

245.7 

237.9 

233.7 

50 

297.8 

285.9 

270.8 

253.2 

240.7 

237.1 

52- 

— 

295.8 

287.3 

272.4 

254.4 

241.6 

237.9 
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Problem 6: Three Dimensions — Figure With Radiative and Flux Surfaces 

A figure with front, back, and two sides has all internal and one external (right side) sur- 
faces radiating with an environment of T e = 255. 6K; the outer surface on the left side is 
exposed to a specified normal flux q = 1576.2 W m' 2 ; and all inner surfaces radiate with 
each other. Tetrahedral elements with grid points as shown have been used to model the 
Figure. Results of the analysis are shown in Table 5. 


SURFACE THICKNESS 0.0254 m 
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Table 5 

Finite Element Results for Problem 6. 


Node 

T at 
3600s 

T at 
7200s 

T at 
18,000s 

T at 
36,000s 

T at 
72,000s 

T for 
Steady State 

1,19. 

310.3K 

334.0K 

404.9K 

509. 8K 

645. 3K 

715. 7K 

2,20 

310.3 

334.1 

404.9 

509.8 

645.3 

715.6 

3,21 

310.4 

334.2 

405.1 

509.9 

645.3 

715.6 

4,16 

318.9 

343.6 

415.0 

521.6 

661.3 

735.1 

5,17 

318.8 

343.5 

414.9 

521.5 

661.2 

735.0 

6,18 

319.1 

343.7 

415.1 

521.7 

661.2 

734.9 

7,13 

330.1 

358.5 

430.2 

537.6 

679.5 

755.1 

8,14 

332.9 

358.4 

430.2 

537.6 

679.4 

754.9 

9,15 

333.0 

358.4 

430.2 

537.5 

679.3 

754.9 

10 

339.6 

365.2 

437.0 

544.4 

686.5 

762.1 

11 

339.5 

365.1 

436.9 

544.4 

686.4 

762.0 

12 

283.6 

364.7 

436.5 

543.9 

685.9 

761.4 

22 

308.7 

332.3 

402.8 

507.0 

640.4 

709.1 

23 

308.7 

332.3 

402.8 

507.0 

640.4 

709.1 

24 

309.1 

332.7 

403.3 

507.4 

640.8 

709.5 

25 

308.7 

332.3 

402.8 

506.9 

640.3 

709.0 

26 

309.1 

332.7 

403.2 

507.3 

640.7 

709.3 

27,45 

309.9 

333.8 

404.6 

509.4 

644.8 

715.1 

28,46 

310.0 

333.8 

404.6 

509.5 

644.8 

715.1 

29,47 

310.1 

333.9 

404.7 

509.6 

•644.8 

715.0 

30,42 

318.9 

343.6 

415.0 

521.6 

661.3 

735.1 

31,43 

318.8 

343.5 

414.9 

521.5 

661.2 

735.0 

32,44 

319.1 

343.7 

415.1 

521.7 

661.2 

734.9 

33,39 

333.4 

358.8 

430.6 

538.1 

680.1 

755.6 . 

34,40 

333.3 

358.8 

430.6 

538.0 

679.9 

755.5 

35,41 

333.4 

358.8 

430.6 

538.0 

679.9 

755.2 

36 

339.7 

365.3 

437.1 

544.5 

686.6 

762.2 

37 

339.6 

365.2 

437.0 

544.5 

686.5 

762.1 

38 

283.7 

364.9 

436.6 

• 544.1 

686.0 

761.6 

48 

308.7 

332.3 

402.8 

506.9 

640.3 

708.9 


\ 
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FINITE ELEMENT COMPUTER PROGRAM 


Input Cards and Data Required 

The following is a description of data necessary for the finite element thermal analysis 
program written for an IBM 360-95 computer. Each data deck must first include card 
types 1 to 5. If radiation boundary conditions exist, these are described by card types 
6 and 7. When time dependent data are used, card type 7 is followed by sets of card types 
8 to 13. (One card type 8, plus combinations of types 9 to 13 must be included for each 
increment where there is a change.) Finally, card types 14 and 15 must be included. 


Card Type 

Data Description 

Data Type 

Field Location 

1 

Title of thermal analysis 

Alphanumeric 

1-80 

Title Card 

problem 

' 


2 

Number of nodes 

Integer (I) 

1-4 

General 

Number of elements 

I 

5-8 

Information 

Number of materials 

I 

9-12 

(one card) 

Dimension of elements 
(2 = triangles, 

3 = tetrahedrons) 

I 

13-16 


Blank spaces 


17-20 


*Number of increments 

I 

21-24 


*Time increment (hr) 

Real (R) 

25-32 


Specified boundary 
code (0 = no, 1 = yes) 

I 

33-36 


Radiation code (0 = no, 

1 = boundary radiation, 

2 = surface radiation) 

I 

37-40 


Convection code (0 = no, 

1 = boundary convection, 

2 = surface convection) 

I 

4144 


Normal flux code (0 = no, 

1 = boundary flux, 

2 = surface flux) 

I 

4548 


Internal heat generation code I 

49-52 


(0 - no, 1 = yes) 




*Note: For steady state results only, place a 0 in field 21-24 (number of increments) and a negative real 
' number in field 25-32 (time increment size). 
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Card Type Data Description Data Type Field Location 


3 

Node number 

I 

1-8 

Nodal 

Information 

x-coordinate (ft) 

R 

9-16 

(one card per 

y-coordinate (ft) 

R 

17-24 

node) 

z-coordinate (ft) 

R 

25-32 


Initial temperature of node 
(°F or 0 R; °R must be used 
if radiative boundary condi- 
tions exist.)* 

R 

33-40 

*Note: If this node temperature is specified (fixed), precede temperature value by a negative sign (-). 

4 

Element number 

I 

1-4 

Element 

*Node i of element 

I 

5-8 

Information 

*Node j of element 

I 

9-12 

(one card per 

*Node k of element 

I 

13-16 

element) 

*Node 1 of element 

I 

17-20 


Material of element 

I 

21-24 


Thickness of element (ft) 

R 

25-32 


Blank spaces 


33-36 


Boundary condition code 
(0 = no boundary conditions, 
1 = boundary conditions) 

I 

37-40 


**Internal heat generation value 

R 

41-48 


(Btu/hr ft 3 ) 




**Convective coefficient h 

R 

49-56 


(Btu/hr ft 2 F) 

**Ambient fluid temperature 

R 

57-64 


(°F or °R) 

tValue of normal flux 

R 

65-72 


(Btu/hr ft 2 ) 



*Note: For triangular elements, describe ij,k in a counterclockwise direction, leaving 1 blank. When a 

boundary of the triangle is subject to boundary conditions, denote the boundary nodes by i and j. 
For tetrahedral elements, again let any surface subject to a boundary condition have nodes i,j,k; 
and number the nodes such that, viewed from point 1, the ij, and k are counterclockwise. 

**Note: These values all are boundary condition variables. Leave the data location blank (or use 0) if the 
information requested is unapplicable to the element described by columns 1-32. 

fNote: The flux in is negative; flux out is positive. 
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Card Type 

Data Description 

Data Type 

Field Location 

5 

Material number 

I 

1-8 

Materia] 

Information 

Conductivity coefficient, k x , 
of material (Btu/hr ft F) 

R 

9-16 

(one card per 
material) 

Conductivity coefficient, ky, 
of material 

R 

17-24 


*Conductivity coefficient K z , 
of material 

R 

25-32 


Heat capacity, pc, of material 

R 

33-40 


*Material internal heat 
generation, Q (Btu/hr ft 3 ) 

R 

41-48 


*Emissivity of material 

R 

49-56 


* Absorptivity of material 

R 

57-64 

*Note: Leave blank if inapplicable to material. If an element 
radiation, 49-64 must be specified. 

or boundary of this material participates in 


If no radiative boundaries exist, omit card types 6 and 7. 


*Temperature of environment R 1-8 

to which continuum is 
radiating 

**Maximum number of iterations I 9-16 

to be performed to attain 
steady state temperatures 

t T , Acceptance criterion for R 17-24 

steady state answers (If left 
blank, = 0.00001 is used.) 


*Note: Leave field 1-8 blank if there is radiation but none with the environment (or with any fixed 
temperature source). 

**Note: This number, N6, should be^> number of increments specified on card type 2, N2, denoting the 
number of increments for which output is desired. However, N6 denotes the maximum number 
of desired iterations (increments) to try to obtain a steady state distribution. Printed results will 
include temperatures at increments 1 through N2 and at steady state (or at N6 if steady state 
has not been attained.) This parameter is required for radiative problems only. 

fNote: Answers are assumed to be at steady state, when for each node 1 1 - T T /T r+c j T | < p. 

This parameter is required for radiative problems only. 


6 

Radiation 
Information 
(one card) 
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Card Type 

Data Description 

Data Type 

Field Location 

7 

* Radiating element 1 

I 

1-4 

Radiation 

* Radiating element 2 

I 

5-8 

Information 

* F l-2 

R 

9-16 

(one card for 

Radiating element 1 

I 

17-20 

five radiating 

Radiating element 3 

I 

21-24 

conditions) 

F l-3 

R 

25-32 


** 




Radiating element n 

I 



Radiating element m 

I 



^ F n-m 

R 



*Note: If 1 and 2 are radiating elements, separate entries must be made for and F^-i > one entr y 
containing “element 1 , element 2, F 12 ” and the other containing “element 2, element 1 , Fj.j 

**Note: Begin new card if needed and proceed as above, using two integer fields of four for identification 
of the two radiating elements, followed by one real field of eight for the form factor value. 

fNote: After all radiative form factors have been read in, place a negative integer in the next integer 
field of four. (I.e., place this negative number where the next element title would be if there 
were another form factor to be described.) This number is a code to denote the end of 
radiative data. 


If no time dependent data exist, omit card types 8—13. 


Increment of first or next I 

change in any input data 
Number of elements with new I 
ambient fluid temperature at 
this increment 

Number of elements with new I 
normal flux value at this 
increment 

Number of elements with new I 
internal heat generation value 
at this increment 

Number of materials with new I 
internal heat generation value 
at this increment 

Number of nodes with new I 

specified temperatures at this 
‘ increment- " - 

*New increment value (if R 

it changes) 


1-4 

5-8 

9-12 

13-16 

17-20 

21-24 

25-32 


*Note: For no change in increment size, leave field 25-32 blank. 
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To describe time changing data for the first (or next) increment, follow card type 8 with 
the appropriate combination of 9, 10, 11, 12, and 13. (I.e., if only the ambient fluid 
temperature is changing at this increment, then only card type 9 is needed to describe it. 
This type, therefore, will follow 8. If, however, there is no change in ambient fluid temper- 
ature, but flux or specified nodal temperature change, card 8 is followed immediately by 
10, then by 13). To describe a change in increment size, only card 8 is needed. 


Card Type 

Data Description 

Data Type 

Field Location 

9 

Number of first element 

I 

1-8 

Time Changing 

subject to new ambient 



Data-Ambient 

fluid temperature 



Fluid Temperature 

New ambient fluid 

R 

9-16 


temperature for element 
Number of second element 

I 

17-24 


subject to new ambient 
fluid temperature 
New ambient fluid 

R 

25-32 


temperature for element 

* 




Number of last element 

I 



subject to new ambient 
fluid temperature 
New ambient fluid 
temperature for element 

R 

1 

i 

1 

*Note: Begin new card if needed and proceed as above, each element having an integer field of 8 for its 
identification and a real field of 8 for its new value. 

10 

Number of first element 

I 

1-8 

Time Changing 

with change in flux value 



Flux Data 

New value of flux 

R 

9-16 


Number of second element 

I 

17-24 


with change in flux value 
New value of flux 

* 

R 

25-32 


Number of last element 

I 



with change in flux value 
New value of flux 

R 


*Note: Begin new card if needed and proceed as before, each element having an integer field of 8 for 
its identification and a real field of 8 for its new value. 
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Card Type Data Description Data Type Field Location 


11 

Number of first element 

I 

1-8 

Time Changing 

with internal heat 



Data-Intemal 

generation change 



Heat Generation 

New value for internal 

R 

9-16 

(Element) 

heat generation 
Number of second element 

I 

17-24 


with internal heat 
generation change 
New value for internal 

R 

25-32 


heat generation 

* 




Same information for last 
element with internal heat 
generation change 



*Note: Begin new card if needed and proceed as before, each element having an integer field of 8 for 
its identification and a real field of 8 for its new value. 

12 

Number of first material 

I 

1-8 

Time Changing 

with internal heat 



Data-Intemal 

generation change 



Heat Generation 

New value of internal 

R 

9-16 

(Material) 

heat generation 
Number of second material 
with internal heat 

I 

17-24 


generation change 
New value of internal 

R 

25-32 


heat generation 

* 




Same information for last 
material with internal heat 
generation change 



*Note : Begin new card if needed and proceed as before, each element having an integer field of 8 for 
its identification and a real field of 8 for its new value. 
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Card Type 

Data Description 

Data Type 

Field Location 

13 

Number of first node 

I 

1-8 

Time Changing 

with new specified 



Data--Specified 

temperature value change 



Nodal Temperature 

New value of nodal 
temperature 

R 

9-16 


Number of second node 

I 

17-24 


with new specified 
temperature 
value change 
New value of nodal 

R 

25-32 


temperature 

* 




Same information for last 
node with specified 
temperature value change 



*Note: Begin new card if needed and proceed as before, each node having an integer field of 8 for its 
identification and a real field of 8 for its new value. 


If data change at a later time increment, return to card type 8. If not, use type 14. 


14 

Integer greater than the 1 

1-4 

End of Time 

maximum increment 


Dependent 

specified on card 2, field 


Data Code 
Card 

21-24 


15 

*Terminator card (i.e., /* for 360-95) 

1-2 

*Note: To submit more than one data deck per run, assemble each deck as above (excluding card type 

15); stack the decks one behind the other; follow the final deck with card 1 5. 
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Computer Program Listing of FEM 


COMMON XIVJZKU2) »X ( 100*3) ,H( 100* 100> *P ( 100* 100) * FORCE U 00) • 

1 AREA (150) »HEAT(150) ♦HCONV ( 1 50 ) *COND (5* 3) *Q (5) *PC (5) *T ( 100 ) ♦ 

2 FO( 100) *F( 100) » AMB( 150) .FLUX (150) « FORM ( 150 » 150) *E(5> .A (5) * 

3 HH (100*100) » TH ( 150 ) * V ( 100) , PER (100) . DET ( A* 4) , VOLUME ( 150) ♦ 

4 STEADY* INDEX (150*4) *NBND1 ( 150 ) ♦ MATERL ( 150 ) . INCR*NEL*NODES*NMAT» 

5 NDIM 

DIMENSION S ( 100* 100) *70(100) *NELEM (5) tMELEM ( 5) *FDR ( 5) * IDENT (20) 
INTEGER*** STEADY 

REAL * 8 DFORCE* DC* DB. DDIF * DDT * DCHK 
C 

NORDER = 100 
10000 CONTINUE 
TIME = 0. 

I NCR = 0 

REAO ( 60 » 99* END=1000 1 ) IDENT 
99 FORMAT (20A4) 

READ (60* 100»END=10001) N00ES»NEL»NMAT*NDIM»MAXINC*DELTA*NC0NST, 

2 NRAD.NCONV.NFLUX.NHEAT 

100 FORMAT (414* I8*F8. 0*514) 

10002 CONTINUE 

WRITE (61* 101) IDENT *NDIM*NEL *NOOES»NMAT *NCONST *NRAD*NCONV*NFLUX* 

1 NHEAT * MAXINC* DELTA 

101 FORMAT ( 1H1/1H0* 45X.20A4/1H0* 

1 / 1H0 *55X ♦ 1 1 HD I MENS I ON =*I3. 

1 /1H0*55X* 20HNUMBER OF ELEMENTS =* 13* 

1 /1H0*55X» 17HNUMBER OF NODES =♦ 13* 

1 /1H0*55X*2 1HNUM8ER OF MATERIALS =* 13. 

1 /1H0.55X.43HC0NSTANT TEMPERATURE INPUT (0=NO. 1=YES) = *11* 

1 / 1H0 » 55X ♦ 68HR ADI AT ION INPUT (0=NO* l = BOUNDARY RADIATION. 2=SURFA 

ICE RADIATION) = *11* 

1 / 1H0 *55X* 71HC0NVECTI0N INPUT (0=NO» l=BOUNDARY CONVECTION. 2=SUR 

IF ACE CONVECTION) = *11* 

1 /1H0.55X*53HFLUX INPUT (0=NO* 1=80UNDARY FLUX. 2=SURFACE FLUX) = 

1 * 11 * 

1 /1H0 *55X»47H INTERNAL HEAT GENERATION INPUT (0=NO* 1=YES) = *11* 

1 / 1 HO ♦ 55X ♦ 35HMAX I MUM NUMBER OF TIME INCREMENTS =*I3* 

1/ 1HO ♦ 55X » 44HT IME INCREMENT (.TRANSIENT* -STEADY STATE) =*F15.5* 

1 / 1HO ) 

C SET SWITCH FOR STEADY STATE (-DELTA) OR TRANSIENT (.DELTA) 

STEADY = 2 

IF (DELTA) 11000*11000.12000 
11000 STEADY = 1 

DELTA = -DELTA 
12000 CONTINUE 
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non oooo o ooo 


NODE INFORMATION 
WRITE (61 * 104) 

104 FORMAT (1H0*50X» 16HN00E INFORMAT ION/ 1H0 * 8X* 1HN, 1 7X, IHX • 14X, 1HY, 

2 14X, 1HZ* 14X. lHT/> 

DO 105 K = 1 * NODES 

READ (60* 102) I , X ( I , 1 ) *X ( I ♦ 2) *X ( I , 3) , T(I) 

102 FORMAT ( I8.4E8*0) 

WRITE (61* 103) I*X(I*1) *X(I*2)*X(I*3)*T(I> 

103 FORMAT (5X*5HN00ES* I3*5X*4F15»5) 

105 CONTINUE 

TABULAR PRINT REGULATOR 
LINE = 60 

LINBLK = (N0DES*5)/5 


ELEMENT INFORMATION 
WRITE (61.113) 

113 FORMAT (1H0/1H0.50X. 19HELEMENT INFORMAT ION/ 1H0 * 5X ♦ 1HE. 5X . IHI » 5X. 

1 1 HJ» 5X ♦ 1 HK* 5X • 1 HL * 3X * 8 HM AT Eft I AL ♦ 7X * 9HTH I CKNESS « 6X * 

1 1H0.8X.7HBC CODE»10X»iHH»llX»lHT» 10X.4HFLUX/) 

DO 112 N = l.NEL 

READ (60 *110) I* ( INDEX ( I * J> ,J*1,4> .MATERL ( I) *TH(I) . 

1 N8ND1 (I) .HEAT (I) *HCONV (I) *AMB(I) » FLUX (I) 

110 FORMAT (614, F8.0* I0.4F6.O) 

WRITE (61*111) I. ( INDEX (I*J)*J=1*4). MATERL ( I ) »TH ( I ) »HEAT ( I ) » 

1 NBND1 (I) » HCONV ( I ) * AMB ( I ) ♦ FLUX (I) 

111 FORMAT (4X. I3*4I6»4X*I4.5X,2F12.3,5X* I5.5X , 3F 12. 3) 

112 CONTINUE 

MATERIAL INFORMATION 
WRITE (61*114) 

114 FORMAT ( 1H0/1HO.50X.20HMATERIAL INFORMATION /1H0* 8X.1HM. 8X*2HKX 

1 * 10X , 2HKY , 10X ♦ 2HKZ » 10X* 2HPC* 10X»1HQ*11X»1HE*11X* 1HA/) 

DO 120 M = l.NMAT 

READ(60* 1 15) I ,COND (1*1)* COND (1*2) *C0ND(I.3> .PC ( I ) , Q ( I ) ,E ( I ) , A ( I ) 

115 FORMAT (I8,7F8#0) 

WRITE (61*116) I * COND (1*1) .COND ( 1 .2) , COND ( I * 3) *PC ( I ) *Q ( I) ,E(I) *A(I) 

116 FORMAT <5X, I5.7F12.3) 

120 CONTINUE 
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n n n non 


radiation Information 

IF(NRAO) 130*130.131 

131 READ (60*1 39) TWALL.MXSSIN .RADTOL 

139 FORMAT (F8.0, I8.F8.0) 

IF <RADTOL .LE. 0.) RADTOL = .00001 
WRITE (61.135) TWALL. MXSSIN ♦ RADTOL 

135 FORMAT (1H0/1H0.50X.1 7HRAD I AT ION FACTORS/ 1HO. 6HTWALL = .F 12. 9. SX. 13HM 

2AX . SS INC . = * 14. 5X. 1 OHTOLERANCE= * F12.6/1H0) 

STEADY = 2 
00 1720 1=1. NODES 
1720 TO ( I ) = ABS (T ( I ) ) 

DO 137 I = l.NEL 
DO 137 J = l.NEL 

137 FORM < I . J) = 0. 

132 READ (60.133) < NELEM ( I ) .MELEM ( I ) . FOR ( I ) . I =1 .5) 

133 FORMAT (5 ( 2I4.F8. 0) ) 

DO 145 1=1.5 
ITEM=I 

IF (NELEM < I ) ) 144,144.136 

136 N=NELEM ( I ) 

M=MELEM( I) 

FORM(N.M) = FOR ( 1 ) 

IF (FORM (N.N) ) 145,138,145 

138 FORM (N.N) = -1. 

145 CONTINUE 
GO TO 146 

144 ITEM=ITEM-1 

IF (ITEM) 130,130.146 

146 WRITE (61. 134) (NELEM ( I ) .MELEM ( I ) .FOR ( I ) .1=1, ITEM) 

134 FORMAT (2X,5(2HF (.12, 1H.,I2,3H) =,F 7.5. 8X1) 

IF (ITEM .EQ. 5) GO TO 132 

130 CONTINUE 

READ (60,140) INC, NAMBS, NFLUXS. NHEATS, NQS. NODS, DELTAN 

140 FORMAT (6I4.F8.0) 

HOW MANY FIXED NODES 

NCONST = 0 
00 150 N = 1 , NODES 
IF (T (N) ) 151,150,150 
151 NCONST = NCONST ♦ 1 
150 CONTINUE 

NVAR = NODES - NCONST 
WRITE (61,154) NVAR, NCONST, NODES 
154 FORMAT ( 1H0, 10X, 5HNV AR = . I3.5X , 7HNC0NST = , 1 3. 5X , 6HN0DES = , 13) 
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CONSTRUCT SIMILARITY TRANSFORMATION MATRIX FOR CONSTANT NOOES 

IF (NCONST .EQ. 0> GO TO 992 
DO 155 I = 1 .NODES 
DO 156 J = 1 ♦ NODES 
156 S<I.J> = 0. 

155 S < I * I ) = 1.0 

NL AST = NODES 
DO 160 N = l.NVAR 
IF < T ( N) > 161*160.160 

161 T (N) = - T ( N) 

IF <T(N) .lE. I.E-7S) T (N) * 0. 

162 IF (T (NLAST ) ) 163.164.164 

163 T < NLAST ) = -T (NLAST) 

IF < T (NLAST) .LE. l.E-75) T (NLAST ) = 0. 

NLAST = NLAST - 1 
GO TO 162 

164 S <N.N) = 0. 

S < NLAST .NLAST ) = 0.. 

S < NLAST *N) = 1. 

S (N.NLAST) = 1. 

NLAST = NLAST - 1 

160 CONTINUE 

NVAR1 = NVAR ♦ 1 
IF (NODES - NVAR 1 ) 992. 167.167 
167 DO 165 N = NVAR1. NODES 
IF<T(N)> 166.165.165 
166 T <N) - -T(N) 

IF (T(N) .LE. l.E-75) T < N ) * 0. 

165 CONTINUE 
C 

992 CONTINUE 
C 

DO 200 I = 1. NODES 
FO(I) = 0. 

DO 200 J = 1. NODES 
H(I.J) = 0. 

P(I.J) = 0. 

200 CONTINUE 

C INITIALIZE THE FORCE VECTOR AND H MATRIX 

NTCD = 0 

15010 DO 199 1=1. NODES 
199 F ( I ) = 0. 
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CALCULATIONS 

DO 2000 N = 1 * NEL 
IF (NTCD »EQ. 2) GO TO 15011 
CALL CONDUC (N*MATERL (N) ) 

C OMIT SPECIFIC HEAT IF STEADY-STATE 

GO T0(15011. 15012) .STEADY 
15012 CALL SPECIF <N, MATERL (N) > 

15011 NN=MATERL(N) 

TOTAL=HEAT (N) ♦ Q(NN) 

IF (TOTAL) 220 » 221 ♦ 220 

220 CALL HE T GEN (N» TOTAL ) 

221 IF (NBNDl(N)) 2000 . 2000 * 224 
224 IF (FLUX (N) ) 226.227*226 

226 CALL FLUXN (NtNFLUX) 

227 IF (HCONV(N)) 228.229.228 

228 CALL CON V EC (NtNCONVtNTCD) 

229 IF (NRAD) 2000.2000.233 

233 IF (FORM (N.N) ) 230*2000.230 

230 DO 232 J = l.NEL 

IF (FORM(N.J)) 231*232*231 

231 CALL RAOIAN(N. J.FORM(N.J) .NN.TWALL.NRAD) 

232 CONTINUE 
2000 CONTINUE 

IF (NTCD .NE. 2) GO TO 234 
NTCD = 0 
GO TO 15113 

234 GO TO (15101* 15103) .STEADY 
15103 DO 201 I=l*NOD£S 

201 FO ( I ) = F ( I ) 

WRITE (61*11119) ( I .VOLUME (I) »I 3 1» NEL ) 

11119 FORMAT (1H0/1H0.50X.26HAREA OR VOLUME OF ELEMENTS 
2 6 (2HA ( » I3.2H) =» F8.3*5X))> 


/(/IX. 
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COMPUTE COEFFICIENT MATRICES <H*2/DELT*P) AND (H-2/DELT*P) 

SKIP IF STEADY-STATE 
15102 FACTOR = 2. /DELTA 
DO 250 I = 1 * NODES 
DO 250 J = I » NODES 
HIJ = H ( I ♦ J) 

PIJ = P(I.J)*F ACTOR 
H ( I * J) = HIJ ♦ PIJ 
H ( J* I ) = H(IiJ) 

P(I*J> = HIJ - PIJ 
250 P(J*I) = P(I*J> 

15101 IF (NCONST .EO. 0) GO TO 329 
DO 321 J=l, NODES 
DO 320 1 = 1, NODES 
V < I ) = 0, 

DO 320 K=l*NODES 

320 V<I> = V < I ) ♦ S(I »K)*H(K*J) 

DO 321 1=1, NODES 

321 H ( I , J) = V ( I » 

DO 326 1=1, NODES 
DO 325 J=l, NODES 
V(J) = 0. 

DO 325 K=l, NODES 

325 V(J) = V(J) ♦ H ( I , K ) * S (K* J) 

DO 326 J=l» NODES 

326 H ( I , J) = V(J) 

329 DO 332 1=1, NODES 

DO 332 J=I, NODES 
HH(IiJ) = HU* J) 

332 HH ( J, I ) = H ( J, I ) 

CALL FACTR (H,PER,NVAR,NORDER* IER) 

IF (IER .NE, 0) WR ITE (6 1 * 1005) IER 

1005 FORMAT (1H0/20H ***** ERROR CODE = ,12, 50H. SEE WRITE UP FOR 

2EXPL ANAT ION. ***** /IHO) 

IF ( NT CO .EQ. 1) GO TO 419 
C 

1 CONTINUE 
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INCR = INCR ♦ 1 
NTCD = 0 

IF (INCR - INO 420.410,420 

410 WRITE (61*500) INCR, NAMBS, NFLUXS, NHEATS*NQS*NODS*DELTAN 
500 FORMAT (1H0/22H CHANGES FOR INCREMENT « 15/lHO* 15* 

1 21H AMBIENT FLUID TEMPS, 5 X, 15, 24H NORMAL FLUX BOUNDARIES* 

2 5X* 15, 24H INTERNAL HEAT ELEMENTS /1H0* 15* 25H INTERNAL HEA 
2T MATERIALS* 5X, 15, 23H SPECIFIED NODAL TEMPS /1H0* 

4F10.5, 20H NEW INCREMENT SIZE /1H0) 

IF (DELTAN) 419* 419* 411 

411 NTCD = 1 

DO 216 1=1 .NODES 
DO 216 U=1 .NODES 
H ( I * J) =0. 

216 P(I,J> = 0. 

DO 217 N= 1 * NEL 

CALL CONDUC (N» MATERL (N) ) 

CALL SPECIF(N*MATERL(N>) 

IF (HCONV(N)) 401.217,401 
401 CALL CONVEC(N.NCONV»NTCD> 

217 CONTINUE 
DELTA = DELTAN 
GO TO 15102 

419 NTCD = 0 

IF (NAMBS) 413*413*412 

412 READ (60*510) (N*AMB(N) * 1=1, NAMBS) 

510 FORMAT (5(18, F8.0)) 

WRITE (61*520) N, AMB (N) 

520 FORMAT ( 1H0 * 4HAMB ( , 1 3* 2H) =* F 12 .5) 

NTCD = 2 

413 IF (NFLUXS) 415*415*414 

414 READ (60*510) (N, FLUX (N) * 1 = 1 .NFLUXS) 

WRITE (61*530) N*FLUX(N) 

530 FORMAT ( 1H0.5HFLUX ( * I3.2H) =*F12.5) 

NTCD = 2 

415 IF (NHEATS) 417*417*416 

416 READ (60*510) (N, HEAT (N) * 1 = 1 .NHEATS) 

WRITE (61*540) N* HEAT (N ) 

540 FORMAT ( 1H0 * 5HHEAT ( * 1 3* 2H) =* F12.5) 

NTCD = 2 

417 IF (NOS) 420*420*418 

410 READ (60.510) (N,Q (N) ♦ 1=1 *NQS) 

WRITE (61*550) N, Q (N) 

550 FORMAT (1H0»2HQ(*I3»2H)=»F12.5) 

NTCD = 2 
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420 CONTINUE 

TIME = TIME ♦ DELTA 
IF (NRAD .NE. 0) NTCD = 2 
IF (NTCO .EO. 2) GO TO 15010 
5113 DO 270 1=1 f NODES 

FORCE (I) = -F(I) -F0(I) 

GO TO <270*15112) .STEADY 
5112 DO 270 J = 1. NODES 

FACTOR = P ( I » J) *T (J) 

FORCE < I ) = FORCE <I) - FACTOR 
270 CONTINUE 

GO TO <152 11*15212) * STEADY 
5212 IF f I NCR - INC) 440.430*440 

430 IF (NODS) 432*432.431 

431 READ (60.510) <N*T (N) ♦ 1=1* NODS) 

WRITE (61*560) N, T (N) 

560 FORMAT ( 1H0 * 2HT ( ♦ I3.2H) =*F 12.5) 

432 READ (60*140) INC.NAMBS*NFLUXS.NHEATS*NQS»NODS* OELTAN 
440 CONTINUE 

COMPENSATE FOR CONSTANT TEMPERATURES 

5211 IF (NCONST ,EO. 0) GO TO 302 

PERFORM SIMILARITY TRANSFORMATION 
DO 327 1 = 1 .NODES 
FO ( I ) = 0. 

V(I> = 0. 

DO 327 J=l. NODES 

FO ( I ) = FO ( I ) ♦ S(I.J) * FORCE < J) 

327 V < I > = V ( I ) ♦ S < I » J) *T ( J) 

DO 328 1=1. NODES 

FORCE (I) = FO ( I ) 

328 T ( I ) = V ( I ) 

K 1 = NVAR ♦ 1 

DO 304 1=1. NVAR 
DFORCE = FORCE (I) 

DC = 0. 

DO 300 J=K1. NODES 
DB = HH(I.J) * T (J) 

300 DC = DC ♦ DB 

DDIF = DFORCE - DC 
304 FORCE ( I ) = DDIF 
302 CONTINUE 

*********************************** 
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TOTAL = 0. 

00 305 I=s 1 fNVAR 

305 TOTAL * TOTAL ♦ (FORCE(I) * FORCE CD) 

TOTAL = SORT (TOTAL) 

IF (TOTAL .GT. I.E-75) 60 TO 310 
TOTAL = 0. 

GO TO 309 

310 DO 306 1 = 1* NV AR 

306 FORCE ( I ) = FORCE < I ) / TOTAL 
EPSI = 2.** (-23) 

CALL RSLMC(HH»H* FORCE* T.NVAR* EPSI »IER*NORDER*V*PER) 

IF (IER .GT. 2) WRITEC61* 1005) IER 
309 DO 307 1=1 *NVAR 

307 T ( I ) = T ( I ) * TOTAL 
DO 308 1=1 *NODES 

308 FOU) = F ( I ) 

*##**##***«****#*»*#**#«•*****#***** 


OUTPUT 

INVERT THE SIMILARITY TRANSFORMATION 
IF (NCONST .EO. 0) GO TO 888 
DO 330 I = 1*N0UES 
V(I) = 0. 

DO 330 J = 1*N0DES 
V < I > = V(I) ♦ S(I*J)*T(J) 

330 CONTINUE 

DO 331 I = 1 * NODES 

331 T ( I ) = V(I) 

OUTPUT 

888 IF (NRAD .EO. 0) GO TO 889 

IF (INCR .LE. MAXINC) GO TO 889 
DO 1700 I=l»NODES 
ODT = T(I) / TO(I) 

DCHK = DABS ( 1 . - DDT) 

CHK = DCHK 

IF (CHK .GT. RADTOL) GO TO 1701 
1700 CONTINUE 

STEADY = 1 
GO TO 889 
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VALUES DID NOT MEET TOLERANCE. THEREFORE PRINTOUT WITH 
ACCEPTABLE TOLERANCE. 

1701 IF ( I NCR - MXSSIN) 1711.1702*1702 

1711 DO 1712 1=1. NODES 

1712 TO ( I ) = TCI) 

GO TO 1 

1702 LINE = LINE ♦ LINBLK ♦ 3 

IF (LINE - 60) 1703.1703.1704 

1704 WRITE (61.997) ID£*IT 
LINE = 10 ♦ LINBLK 

1703 TOL = CHK 

DO 1706 J=I .NODES 
DDT = T (J) / TO ( J) 

DCHK = DABS (1. -DDT) 

CHK = DCHK 

IF (CHK .GT. TOL) TOL = CHK 
1706 CONTINUE 

WRITE (61.1705) TOL 

1705 FORMAT ( 1H0/3X * 78HSTEADY STATE TOLERANCE NOT MET. ..AT MAXIMUM SS I 

CNCREMENT. ABS ( 1-T (N) /T (N-l) ) =, E12.6. 1H./1H0) 

STEADY = 1 
GO TO 15122 

889 LINE = LINE ♦ LINBLK ♦ 3 
IF (LINE - 60) 901.901.900 

900 WRITE (61.997) IDENT 

997 FORMAT ( 1H1/1H0.25X.20A4/1HO) 

LINE =10 ♦ LINBLK 

901 GO TO ( 15121. 15122) .STEADY 

15121 WRITE (61.998) INCR 

998 FORMAT ( 1H0/30X . 45HSTE ADY-ST ATE NODAL TEMPERATURES INCREMENT #. 

1I4/1H0) 

GO TO 15123 

15122 WRITE (61.999) INCR. TIME 

999 FORMAT (10H0INCREMENT. 15. 5X.9HAT TIME =.F 12. S. 5X. 50 ( 1H-) /) 

15123 WRITE (61.1000) ( I . T ( I ) . I = 1 .NODES ) 

1000 FORMAT (5(6X»2HT (.I2.3H) =»F12.4> ) 

GO TO ( 10000.3) .STEADY 
3 IF ( INCR-MAXINC) 1.2.2 
C CHANGE TRANSIENT TO STEADY-STATE 

2 IF (NRAO .LE. 0) GO TO 1710 
DO 1709 1=1. NODES 

1709 TO ( I ) = T ( I ) 

GO TO 1 

1710 STEADY = 1 
GO TO 992 

10001 STOP 
END 
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SUBROUTINE CONDUC<N.M) 

COMMON XI.XJ.XK*XL.YI.YJ.YK«YL.ZI.ZJ*ZK.ZL.X<100.3)«H<i00.100). 

1 P(IOO.IOO) .FORCEUOO) ♦ 

1 AREA ( ISO) .HEAT <150) .HCONV< 150) *C0ND<5.3> »Q (5) .PC <5> .T < 100) ♦ 

2 FO < 100 ) . F ( 100) * AMB ( 1 50 ) .FLUX 1 1 50) .FORM <150.1 50 ).E(5)*A(5). 

3 HH <100.100) *TH<150) . V < 100) .PER <100) .0ET<4.4) .VOLUME < 150) . 

4 STEADY. INDEX <150.4) . NBND1 < 150) ♦ MATERL < 150 ) . INCR.NEL. NODES. NMAT. 

5 NDIM 

INTEGER * 4 STEADY 

CONDUCTIVITY MATRIX 

I = INDEX <N» 1 ) 

J = INDEX (N» 2) 

K = INDEX <N.3) 

XI = X<I,1> 

XJ = X(J.l) 

XK = X <K • 1 ) 

Y I = X ( I . 2) 

YJ = X ( J. 2) 

YK = X <K *2) 

IF <NDIM ,EQ. 3) GO TO 100 

AREA <N) = A8S <X J*YK - XK*YJ - XI*<YK-YJ> ♦ YI* <XK-XJ) ) /2. 

VOLUME (N> = AREA <N) 

BI = YJ-YK 
BJ = YK- Y I 
BK = YI-YJ 
Cl = XK-XJ 
CJ = XI-XK 
CK = XJ-XI 

CX = C0ND<M,1)*TH<N)/<4.*AREA<N) ) 

CY = C0ND<M.2)*TH<N)/(4.*AREA<N) ) 


H< 

1. 1) 

3 

H < I ♦ I ) 

♦ 

BI*CX*BI 

♦ 

CI*CY*CI 

H( 

I. J) 

X 

H < I . J) 

♦ 

BJ*CX*BI 

♦ 

CJ*CY#CI 

H < 

I.K) 

= 

H < I » K ) 

♦ 

BK*CX«BI 

♦ 

CK*CY*CI 

H< 

J. J) 

= 

H< J. J) 

♦ 

BJ*CX*BJ 

♦ 

CJ*CY*CJ 

H < 

J.K) 


H< J.K) 

♦ 

8K*CX*B J 

♦ 

CK*CY*CJ 

H ( 

K.K) 

= 

H<K.K) 

♦ 

BK*CX*BK 

♦ 

CK*CY*CK 

H < 

J.I) 

= 

H< I. J) 





H < 

K.I) 

= 

H < I . K ) 





H ( 

K. J) 


H < J ♦ K ) 






RETURN 
100 CONTINUE 
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CONDUCTIVITY MATRIX FOR THE 3 DIMENSIONAL TETRAHEDRON 

L=INDEX (N*4> 

XL = X (L » 1 ) 

YL = X (L»2) 

Z I = X(It3) 

ZJ = X ( J» 3) 

ZK = X (K»3) 

ZL = X (L* 3) 

DET(ltl) = XJ*(YK*ZL-YL*ZK) - YJ* (XK*ZL“XL*ZK) ♦ ZJ* < XK*YL-XL*YK) 
DET (2*1) = -YK*ZL ♦ YL*ZK ♦ YJ*(ZL - ZK) - ZJ«<YL - YK) 

DET (2» 2) =• ( -YL*ZI ♦ YI*ZL ♦ YK*(ZI-ZL> - ZKMYI-TL)) 

DET (2*3) = -YI*ZJ ♦ Y J*Z I ♦ YL*(ZJ - ZI) - ZL*(YJ - YI) 
DET(2*4)=-(-YJ*ZK ♦ YK*ZJ ♦ YIMZK - ZJ) - ZI«(YK * YJ) ) 

DET (3» 1 ) =- ( XJ* ( ZL - ZK) -<XK*ZL - XL*ZK) ♦ ZJ*<XK - XL)) 

DET (3t 2) = XK* (ZI - ZL) - (XL*ZI - XI*ZL) ♦ ZK*(XL - XI) 

DET (3*3)=-(XL*(ZJ - ZI) - (XI*ZJ - XJ*ZI) ♦ ZL*<XI - XJ)) 

DET(3»4)= XI* (ZK - ZJ) - <XJ*ZK - XK*ZJ) ♦ZI*<XJ - XK) 

DET (4» 1 ) = -XJ* ( YK— YL) ♦ YJMXK - XL) - XK*YL ♦ XL*YK 
DET(4,2)=-(-XK*(YL - YI) ♦ YK*<XL - XI) - XL*YI ♦ XI*YL) 

DET (4» 3) = -XL* ( Y I -YJ) ♦ YL*(XI - XJ) - XI*YJ ♦ XJ*YI 
D£T (4»4> =-<-XI* (YJ - YK) ♦ Yl*(XJ - XK) - XJ*YK ♦ XK*YJ ) 

IF (NBNDl(N) .EQ. 0) GO TO 200 

AREA (N) = SORT ( . 25* (DET ( 2»4) **2 ♦ D£T(3,4)**2 ♦ DET (4,4) **2) ) 

200 VOLUME (N) = ABS ( (DET ( 1 t 1 ) ♦DET (2t 1 > *XI*DET < 3» l > *YI *DET <4* 1 ) *Z I > /6. ) 
CW = 36.*V0LUME(N) 

CX = COND (M» 1 ) /CW 
CY = C0ND(M»2)/CW 
CZ = COND (M» 3) /Crt 
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H(I.l) * H ( 1 , 1 ) ♦ OET 12* l)«*CX*OET (2*1) ♦ OET (3* 1>*CY*0ET (3*1) 

2 ♦ DET (4, 1) *CZ*OET (4, 1) 

H(ItJ) = H ( I » J) ♦ DET (2*2) *CX*DET (2,1) ♦ DET (3,2> *CY*DET <3. 1 ) 

2 ♦ DET (4,2) *CZ*DET (4, 1) 

H(I.K) = H ( I »K) ♦ DET ( 2. 3) *CX*OET (2*1) ♦ DET (3.3) *CY*DET <3, l > 

2 ♦ DET (4,3) *CZ*DET (4, 1) 

H ( I »L ) = H ( I ,L ) ♦ DET (2,1) *CX*DET(2,4) ♦ OET (3» 1 > *CY*OET < 3,4) 

2 ♦ DET(4,1)*CZ«DET<4,4) 

H ( J, J) = H(J.J> ♦ DET (2»2>*CX*DET (2,2) ♦ DET <3.2> *CY*DET (3,2) 

2 ♦ DET (4,2) *CZ*OET (4,2) 

H ( J, K ) = H(J,K) ♦ DET (2.3)*CX*0ET (2,2) ♦ DET (3.3) *CY*0£T (3.2) 

2 ♦ DET (4,3) *CZ»DET (4,2) 

H(J,L) = H ( J, L ) ♦ DET (2,2)*CX*DET(2,4) ♦ DET (3,2) *CY*DET (3,4) 

2 ♦ DET (4,2)*CZ*DET (4,4) 

H (K,K) = H(K,K) ♦ DET(2»3)*CX*DET<2,3) ♦ DET (3,3) »CY*DET (3,3) 

2 ♦ DET (4,3) *CZ*DET (4,3) 

H(K,L) = H(K,L) ♦ DET ( 2, 3) *CX*DET ( 2, 4) ♦ OET (3,3) *CY*OET (3.4) 

2 ♦ DET (4,3)*CZ»DET (4,4) 

H<L,L) = H(L.L) ♦ 0ET(2,4)*CX*DET<2,4) ♦ OET (3,4) *CY*OET (3.4) 

2 ♦ DET (4.4) *<~ 7 *DET (4.4) 

H ( J. I ) = H(I.J) 

H ( K ♦ I ) = H(I.K) 

H(L,l) = H ( I , L) 

H ( K ♦ J) = H ( J » K ) 

H(l_,J) = H(J,L) 

H(L,K) = H(K.L) 

RETURN 

END 
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SUBROUTINE SPECIF(N.M) 

SPECIFIC HEAT MATRIX 

COMMON XX (4) »Y <4) .Z<4> .X(100.3) «H(100» 100) *P<100» 100) .FORCE (100) ♦ 

1 AREA < 150) » HEAT (150) * HCONV (150) .C0ND<5.3) .Q <5> .PC ( 5) • T ( 100 ) » 

2 FO(IOO) »F < 100) »AMB(150) .FLUX ( 150) .FORM ( 150 » 150 ) »E(5)«A(5) * 

3 HH <100. 100) »TH(150> ♦ V < 100 ) . PER < 100 > .DETC4.4) .VOLUME (150) . 

4 STEADY. I NOEX< 150.4) .NBNDl (150) .MATERL < 150) » INCR.NEL* NODES .NMAT. 

5 NDIM 

I = INDEX (N. 1 ) 

J = INDEX (N» 2) 

K = INDEX <N.3) 


IF (NDIM 

.EQ. 3) 


GO TO 100 

FACTOR 


PC < M) * AREA (N)°TH(N) / 12 

F ACT 02 

= 

FACTOR 

♦ 

FACTOR 

10 P(I*I) 


P(I.I) 

♦ 

F ACT 02 

P ( I • J) 

= 

P < I « J) 

♦ 

FACTOR 

P ( I » K) 

= 

P( I.K) 

♦ 

FACTOR 

P ( J* J) 


P< J. J) 

♦ 

FACT02 

PC J.K) 

s 

P ( J.K) 

♦ 

FACTOR 

P(K.K) 


P(K.K) 

♦ 

FACT02 

P(J.I) 


P(I.J) 



P(K.I) 


P(I.K) 



P(K.J) 

RETURN 


P ( J.K) 




100 CONTINUE 

L=INDEX (N.4) 

FACTOR = PC ( M ) * VOLUME ( N) * # 05 
FACT02 = FACTOR ♦ FACTOR 
P(I.L) = P(I.L) ♦ FACTOR 

P(J.L) = P(J.L) ♦ FACTOR 

P(K.L) = P(K.L) ♦ FACTOR 

P(L.L) = P(L.L) ♦ FACT02 

P(L.I) = P < I .L) 

P(L.J) = P < J.L) 

P(L.K) = P <K.L) 

GO TO 10 
END 
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SUBROUTINE RADIAN (N*M*FNM»MAT *TWALL*NRAD) 

C 

CRADIATION NORMAL TO BOUNDARY 
C 

COMMON X I Y JZK ( 12)*X(100.3)*H(100.100).P(100.100) ♦ FORCE (100) » 

1 AREA (150) » HEAT (150) .HCONV( 150) .C0ND(5.3) *0(5) .PC < 5) »T (100) ♦ 

2 FO(10Q)*F(100> * AMB( 150) .FLUX ( 150 >» FORM < 150 * 150 ) *E(5) .A(5) . 

3 HH< 100 * 100) » TH ( 150) ♦ V ( 100) .PER ( 1 00 ) ♦ DET (4*4) . VOLUME ( 150 ) ♦ 

4 STEADY* INDEX (150*4) * NBND1 ( 150) ,MATERL( 150) * INCR*NEL*NODES»NMAT . 

5 NDIM 

REAL * 8 DF1.DF2.DF3.DT4 
I = INDEX (N* 1 ) 

J = INDEX (N* 2) 

IF (NRAO .EO. 2) 60 TO 100 

XL 1 a SORT ( (X(I»1)-X(J*1> )**2 ♦ (X(I»2)-X(J»2) )**2) * TH (N) 

IF (N.NE.M) GO TO 20 

T41 a ( .005* (TCI) ♦ T(J> )>**<* * .0857 * XL 1 * E (MAT ) 

IF (FNM) 5*6.6 

5 NOWALL = 0 
GO TO 10 

6 NOWALL = 1 
10 DF1 a F ( I) 

DF 2 - F ( J) 

DT 4 a T41 
DF 1 a DF 1 ♦ DT 4 
DF2 a DF2 ♦ DT 4 
F ( I ) = DF 1 
F ( J) = DF 2 

IF (NOWALL .EO. 0) RETURN 
T41 a (TWALL * .01) ** 4 
GO TO 30 
20 K= INDEX ( M . 1 ) 

L=INDEX (M.2) 

T 1 = .005 * (TOO ♦ T (L> ) 

T41 = T1 ** 

30 T41 a -T41 * .0857 * XL1 * A ( MAT ) * FNM 
NOWALL = 0 
GO TO 10 
100 CONTINUE 

K= INDEX (N» 3 ) 

A 1 = AREA(N) 

T 1 = (T ( I ) ♦ T ( J ) ♦ TOO) / 3. 

IF <N .NE. M) GO TO 220 
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o n 


IF N EQUALS M» CALCULATE RADIATION EMITTED ONLY PLUS THAT ABSORBED 
CFROM KNOWN SOURCE. 

C 

T41 = ( (.01*T1>**4)* .1714 « A1 * E (MAT ) / 3. 

IF (FNM) 205*206*206 

205 NOWALL = 0 
GO TO 210 

206 NOWALL = 1 
210 DF1 = F ( I ) 

OF 2 = F ( J) 

OF 3 = F (K) 

DT4 = T41 
OF 1 = DPI ♦ DT 4 
DFa = DF2 ♦ DT 4 
DF3 = DF3 ♦ DT4 
F C I » = DF1 
F ( J) = DF2 
F (K) = DF3 

IF (NOWALL .EQ. 0) RETURN 
T41 = (TWALL * .01) ** 4 
GO TO 230 

220 II = I NDEX (M* 1 ) 

JJ = INDEX (M.2) 

KK = INDEX ( M* 3) 

T 1 * .01 * ( T (II) ♦ T ( JJ) ♦ T (KK ) ) / 3. 

T41 = T1 ** 4 

230 T41 = -T41 * .1714 * A1 * A (MAT) * FNM / 3. 

NOWALL = 0 
GO TO 210 
END 
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SUBROUTINE HETGEN (N.TOT AL> 

C 

C INTERNAL HEAT GENERATION - FORCE 

C 

COMMON XIYJZKU2) ♦ X ( 1 00.3) «H ( 100. 100) »P ( 100 ♦ 100 ) .FORCE ( 100) ♦ 

1 AREA (150) .HEAT (150) .HCONV(ISO) »C0ND(5»3) .Q (5) .PC < 5) *T 1 100) * 

2 FO(IOO) »F (100) ♦ AMB ( 150 > .FLUX ( 150 ) .FORM < 150 « 150) .E<5) ,A(5) . 

3 HH(IOO.IOO) *TH(150) ,V(100) .PER (100) »DET (A, 4) .VOLUME! 150) . 

4 STEADY . INDEX ( 150.4) .N6ND1 ( 150) .MATERL ( 150) ♦ INCR.NEL. NODES. NMAT. 

5 NDIM 

I = INDEX(N.I) 

J = INDEX (N.2) 

K * INDEX (N» 3) 

IF (NDIM .EQ. 3) GO TO 100 

FACTOR = -TOTAL * AREA (N) * THIN) /3. 

10 F ( I ) = F ( I ) ♦ FACTOR 
F(J) = F < J) ♦ FACTOR 
F (K) = F (K) ♦ FACTOR 
RETURN 
100 CONTINUE 

L = INDEX (N.4) 

FACTOR = -T0TAL*V0LUME(N)*.25 
F <L) = F <L) ♦ FACTOR 
GO TO 10 
END 
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SUBROUTINE FLUXN (N»NFLUX) 


FLUX NORMAL TO SURFACE - FORCE 

COMMON XIYJZK(12) *X ( 100. 3) *H ( 100* 100) t P ( 1 00 » 10 0 ) .FORCE < 100) ♦ 

1 AREA (ISO) »H£ATU50> »HCONV ( 150) tC0ND(5*3> *Q (5) . PC (5) ♦ T ( 100 ) t 

2 F0(100) »F (100) »AMB( ISO) t FLUX ( 150 ) *FORM ( 150 » 150 ) «E(5> *A(5> ♦ 

3 HH(100»100) *TH(150)*V(100) t PER (100) , DET ( k% A) ♦ VOLUME ( ISO > * 

A STEADY * INDEX (150*4) tNBNDl (150) » MATERL ( 150) * INCR»NEL* NODES tNMAT* 
5 NDIM 

NR = INDEX (N$ 1) 

NS = INDEX (N»2) 

IF (NFLUX .EQ. 2) GO TO 200 

D= SORT ( (X (NR* 1 ) -X (NS* 1 ) )**2 ♦ (X (NRt2)-X (NS« 2) ) **2) * TH(N) 
FACTOR = FLUX (N) ° D * .5 
F (NR) = F (NR) ♦ FACTOR 
F (NS) = F (NS) ♦ FACTOR 
RETURN 

200 NT = INDEX (N* 3) 

FACTOR = (FLUX (N) * AREA (N) ) / 3. 

F (NR) = F (NR) ♦ FACTOR 
F(NS) = F ( NS) ♦ FACTOR 
F ( NT ) = F ( NT) ♦ FACTOR 
RETURN 
END 
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SUBROUTINE CONVEC (N.NCONV*NTCD> 


CONVECTION ACROSS BOUNDARY ♦ AMBIENT 

COMMON XIYJZM 12) *X (100*3) *H(100*100) .P( 100* 100) *FORCE(100> ♦ 

1 AREA (150) * HEAT (150) .HCONV ( 150) .COND (5.3) . Q ( 5) .PC < 5) . T ( 1 00 ) . 

2 FO(IOO) »F(100) .AMB( 150). FLUX (150) .FORM (150 *150) *E(5).A(5). 

3 HH(100» 100) *TH( 150) »V(100) ♦ PER (100) .GET (4, A) .VOLUME (150) ♦ 

4 STEADY. INDEX (150*4) .NBNDl (150) .MATEPL(lSO) * INCR.NEL.NODES*NMAT* 

5 NDIM 
I = INDEX (N. 1 ) 

J = INDEX (N« 2) 

IF ( NCONV .EQ. 2) GO TO 100 

HI = SORT ( (X ( I « 1 ) - X(J.1>>**2 ♦ (X ( I . 2) - X ( J, 2) ) **2> *TH (N) 

2 *HC0NV(N)/6. 

IF (NTCD .EQ. 2) GO TO 2 
H2 - 2 » ° HI 
H(I,I) = H ( I » I ) ♦ H2 
H(J.J) = H(J.J) ♦ H2 
H ( I , J) = H ( I * J) ♦ HI 
H ( J, I ) = H ( I * J) 

IF (NTCD .EQ. 1) RETURN 
2 H4 = -H1*AMB (N) *3. 

F ( I) = F ( I) ♦ H4 
F ( J) = F ( J) ♦ H4 
RETURN 
100 CONTINUE 
C 

CONVECTION OVER SURFACE 
C 

K = INDEX (N. 3) 


HI 

= HCONV 

(N) * 

AREA (N) 

IF 

(NTCD 

• 

EQ. 2) 

GO TO 

H2 

= HI ' 


2. 


H ( I 

.1) = 

H 

(I. I) 

♦ H2 

H ( J 

♦ J) = 

H 

< J. J) 

♦ H2 

H (K 

. K ) = 

H 

(K.K) 

♦ H2 

Hd 

*J) = 

H 

(I. J) 

♦ HI 

H ( I 

.K> = 

H 

( I ♦ K ) 

♦ HI 

H( J 

♦ K) = 

H 

( J.K) 

♦ HI 

H ( K 

. J) = 

H 

( J.K) 


HOC 

*1) = 

H 

(I.K) 


H ( J 

.1) = 

H 

(I. J) 


IF 

(NTCD 

• 

EQ. 1) 

RETURN 


200 H4 = HI * 4. * AMB(N) 
F ( I ) = F(I) - H4 

F ( J) = F ( J) -H4 
F (K) = F (K) - H4 
RETURN 
END 
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